Prediction criterion and numerical validation for the interaction between hydraulic fractures and bedding planes

Shale is a kind of sedimentary rock with an obvious bedding structure. The effect of the bedding plane on hydraulic fracture initiation, propagation, and complex fracture network formation is remarkable and a major problem in hydraulic fracturing and shale oil and gas development. In this study, a criterion is established to predict the evolution behavior of hydraulic fractures (HF) under different confining pressure differences and intersection angles. This criterion is intended to predict the types of interaction between HFs and bedding planes (BPs): penetrating, slipping, or dilating. The dependence of crossing on the intersection angle and the principal stress difference is quantitatively presented using the criterion. Meanwhile, 20 simulations with principal stress differences of 2, 4, 6, and 8 MPa and intersection angles of 30°, 45°, 60°, 75°, and 90° were simulated using the RFPA2D-Flow code. Simulation results exhibit good agreement with the criterion results for a wide range of angles. The investigation showed that HFs tend to penetrate BPs under high confining pressure differences and intersection angles and open BPs under low confining pressure differences and intersection angles. In addition to the above two forms, HFs slip due to shear. The criterion can provide relevant reference about the formation of complex fracture networks in shale layers.


Introduction
Shale oil and gas are a kind of clean, efficient, and environmentally energy.As one of the unconventional oil and gas resources, shale is receiving increasing attention around the world because it is one of the most realistic alternative resource replacement for conventional oil and gas [1][2][3][4][5].As a kind of sedimentary rock, shale is characterized by significant bedding and its formations are known to have a fine layered structure [6].The existence of bedding planes (BP) affects the stress distribution, thus affecting the fracture propagation trajectory [7].Many studies have shown that the existence of bedding promotes the formation of fracture networks in shale reservoirs [8][9][10][11][12][13], thus promoting shale reservoir permeability and oil and gas production [14][15][16].Since 1947, the year of the first experimental hydraulic fracturing treatment in the United States in the Hugoton gas field in Grant County, Kansas [17], hydraulic fracturing treatments have become an essential part in gas and oilfield development [18][19][20].Therefore, the investigation of the interaction between hydraulic fractures (HFs) and BPs is of great significance for the complex network formation [13], the reconstruction of shale reservoir, and even the exploitation of shale oil and gas.
In recent years, some scholars studied hydraulic fracturing problems, especially the HF evolution in discontinuous shale reservoirs and have drawn some valuable conclusions.Weak discontinuities, such as natural fracture(NF) and moderate development bedding, are beneficial to oil-gas recovery because the interaction between HFs and discontinuous is an important factor of the hydraulic fracture network's complexity [12,[21][22][23].The interaction between HFs and discontinuities(i.e.BPs in this study) can be observed in three main modes [6,7,11,19,24,25]: 1. HF opens the BP and propagates along the BP (dilate); 2. HF is captured by the BP, and then the BP shear slips (slip); 3. HF passes through the BP (penetrate), as shown in Fig 1 .Through numerous theoretical analyses, physical experiments, and numerical simulations, several criteria have been proposed to predict the occurrence of these three interaction modes, such as the Blanton [26], Warpinski [27], Renshaw [28], Zhou [25], Gu [19], Cheng [29], Sarmadivaleh [30], and Gong [31] criteria.Some of these criteria depend on the output parameters in the construction process, some only consider whether HF penetrate the discontinuities, some are extremely complex to use directly, and most studies focus on the interaction between HF and NF.
This paper proposes a usable criterion of the interaction between HF and BP on the basis of fracture mechanics to predict the dilating, slipping, or penetrating behavior of HF after encountering BP using some parameters as input.The HF evolution of the shale layer is simulated at different principal stress differences and intersection angles to verify the criterion.

Interaction criterion between HF and BP
The interaction between a HF and a BP is a complex process.On the basis of the experiment results and theoretical analysis, three modes may form during the propagation of the HF toward the BP and beyond.First, when the fracture tip reaches the interface, the BP will dilate (Fig 1A ) and the HF will propagate along the BP if the pressure in the HF exceeds the sum of the normal stress acting on the BP and the tensile strength of the bedding.On the contrary, if the pressure in the HF is less than the sum of the normal stress acting on the BP and the tensile strength of the bedding, two possible outcomes form the interaction between the HF and the BP, namely, slipping (Fig 1B)

and penetration (Fig 1C).
To derive the prediction criterion, a rectangle model is set up.As shown in Fig 2, the minimum principal stress (σ 3 ) is applied in the vertical direction, and the maximum principal stress (σ 1 ) is applied in the horizontal direction.The BP approaches HF at an arbitrary angle (θ).The In Region Ⅰ, The stress in the HF and the BP is considered uniform.Therefore, Then, solutions can be obtained by solving the above equilibrium Eqs (1 and 2): When the fracture tip reaches the interface, the BP will dilate if the pressure in the HF exceeds the sum of the normal stress acting on the BP and the tensile strength of the bedding.The dilation condition is p>σ n1 +T O or p>σ n2 +T O .Given that σ n1 >σ n2 , the dilation condition is simplified to where T O is the tensile strength of the bedding.Then, Eq 5 is substituted into Formula 6 to obtain the following: When the fracture tip reaches the interface, shear slip may occur on the BP if the pressure in the HF is less than the sum of the normal stress acting on the BP and the tensile strength of the bedding.The slip conditions are p<σ n1 +T O and Given that σ n1 >σ n2 and τ 2 >τ 1 , the slip condition is simplified to where S O is the shear strength of the bedding, and f is the friction coefficient of the bedding.Then, Eq 5 is substituted into Formulas 8 and 9 to obtain the following: When the fracture tip reaches the interface, the HF may cross the BP if the pressure in the HF is less than the sum of the normal stress acting on the BP and the tensile strength of the bedding and shear slip does not occur.The penetration conditions are p<σ n1 +T O and p<σ n2 +T O , τ 1 <S O +f(σ n1 −p) and τ 2 <S O +f(σ n2 −p).Given that σ n1 >σ n2 and τ 2 >τ 1 , the penetration condition is simplified to Formula 8 and 12 By substituting Eq 5 into Formulas 8 and 12, Formula 10 and be obtained.
In fracture mechanics, the propagation of HF is controlled by the fracture toughness (K IC ) .The half length of the HF is a, and stress intensity factor K I ¼ s ffi ffi ffi ffi ffiffi pa p . When the HF approaches the BP, is always satisfied.By substituting Eq 14 into Formulas 7, 10, 11, and 13, the interaction criterion can be obtained.
Dilation criterion Slip criterion Penetration criterion The criterion is the function of the principal stress difference and the intersection angle.Parameters a,b,To,f, and K IC in the criterion can be input according to the actual situation.

Model set
The numerical simulation of the interaction between the HF and the BP can be divided into three parts.The first part refers to the experiment in [32] for verifying the reliability of the RFPA-Flow code [9,10,[33][34][35].The second part refers to the verification of the above prediction criterion.The third part refers to the investigation on the influence of the absolute value of principal stress on the interaction.As shown in Fig 3, the dimension of the model is 300 mm × 300 mm, which is divided to form a 300 × 300 mesh.The spacing of the parallel bedding is 37.5 mm, and the thickness of the bedding is 1 mm (single element).The wellbore is located in the center of the model, with a radius of 10 mm.In the first part, the intersection angle (θ) is 90˚and the principal stress difference is 0.5 MPa, which are according to the experiment.In the second part, the intersection angle is 30˚, 45˚, 60˚, 75˚, and 90˚at principal stress differences of 2(σ 1 = 4MPa, σ 3 = 2MPa), 4(σ 1 = 6MPa, σ 3 = 2MPa), 6(σ 1 = 8MPa, σ 3 = 2MPa), and 8 MPa(σ 1 = 10MPa, σ 3 = 2MPa).In the third part, the intersection angle and the principal stress difference are same (θ = 30˚, 4σ = 2MPa,) but the maximum and the minimum principal stress are different (σ 1 = 4MPa, σ 3 = 2MPa; σ 1 = 6MPa, σ 3 = 4MPa; σ 1 = 8MPa, σ 3 = 6MPa; σ 1 = 10MPa, σ 3 = 8MPa; σ 1 = 12MPa, σ 3 = 10MPa; σ 1 = 14MPa, σ 3 = 12MPa; σ 1 = 16MPa, σ 3 = 14MPa; σ 1 = 18MPa, σ 3 = 16MPa).The HFs involved in this study initiated and propagated in the direction of the maximum principal stress, so the intersection angle (θ) is the angle between the BP and the direction of the maximum principal stress.A plane-strain component is employed for calculation.
Rock mass is a highly heterogeneous material containing geological discontinuities such as faults, bedding planes, joints and microcracks [6].Therefore, the consideration of heterogeneity is particularly important in the numerical simulation of rock.The mesoscopic elements in RFPA models are assumed to be isotropic and homogeneous, and their mechanical properties (e.g., Young's modulus, Compressive strength, among others) are also assumed to be linear.The mesoscopic elements are statistically distributed (e.g., normal, Poisson, and Weibull distributions) to describe the mechanical properties of the nonlinear macro models.In the statistical distribution density of the mechanical parameters of the mesoscopic elements in RFPA codes, a heterogeneity index m is used to describe the heterogeneity of solid materials.A high m value indicates the presence of highly homogeneous materials whereas a small m value denotes the existence of highly inhomogeneous materials.The meso-mechanical parameter values which are input in modeling should be obtained through the simulations of mechanical property experiments by RFPA 2D-Basic.During this period, the suitable mesoscopic values would be got when macro mechanical properties of the simulations are close to the experimental values through multiple adjustments of the mesoscopic mechanical parameters.

Numerical simulation of layered shale fracturing experiment
According to the experiment by [32], the value of σ 1 (corresponding to σ V in the experiment) is 20 MPa and that of σ 3 (corresponding to σ H in the experiment) is 19.5 MPa.In the experiment, the direction of perforation is normal to the bedding planes.The vertical in situ stress (σ V ) is loaded normal to the bedding planes and the greatest horizontal in situ stress (σ H ) is loaded along the beddings.The list of rock and bedding parameters used in the simulations is presented in Table 1.The simulation conditions are summarized in Table 2.As shown in Fig 5, the HF initiated and propagated in the vertical direction due to the influence of the insitu stress.Therefore, the principal stress difference is 0.5 MPa, and the intersection angle θ is 90˚.
Acoustic emission (AE) monitoring and the specimen cut were used in the hydraulic fracturing experiment to analyze the results in detail.The AE energy is low when the HF propagates along weak beddings; thus, the local complex sub-fractures are difficult to capture by AE monitoring [32].RFPA simulation solves this problem well because it can simulate the failure process of heterogeneous and permeable geo-materials and capture the AE events throughout the entire failure process, and it even determines whether the AE event belongs to tensile or compressive failure.
The spatio-temporal evolution of HFs is expressed by the pore pressure field, the minimum principal stress field, and the AE events distribution in Fig 4 .As shown in the figure, the HFs initiated at the top and bottom of the wellbore and propagated at an angle of approximately 20ẘ ith the greatest in situ stress, mainly because of the low confining pressure difference.Then, the HFs approached the BP at an intersection angle of 90˚for a period of growth; the top HF penetrated one BP and then the next BP shear slip, while the bottom one slipped directly.Thus, slipping was mainly caused by tensile-shear failure induced by the presence of BPs, which is weaker than rock the matrix.The HF at the bottom had a slow propagation speed than that at the bottom due to shear slip, so the fracture was shorter.Furthermore, the HFs that extended at low speed communicated with the weak beddings by slipping and penetration, thus producing a complex fracture network.The green region around the fracture tip in the minimum principal stress pictures show that the HFs propagation is due to the tensile stress.However, the white circle in the BP in the AE events distributions (a3 and c3 in Fig 4), which indicates that the Table 1.Physicomechanical parameters of the rock and bedding in the simulations.

Unit Rock matrix Bedding
Heterogeneity index (m) 3 3 Young's modulus (E) GPa 30(meso value is 37) 3(meso value is 3. compression failure shows that tensile failure (red circle in AE events distributions) is the main but not the only form of failure in hydraulic fracturing in the shale layer.
In the simulation, plane strain calculation was used to reduce the influence of the twodimensional model.The intersection between the HF and the BP, such as penetration and the slip in the simulation during the HF evolution, is basically consistent with the experiment (Fig 5).The result also verifies the reliability of the RFPA code.

Numerical simulations for verifying the prediction criterion
In this part, twenty simulations with principal stress differences of 2, 4, 6, and 8 MPa and intersection angles of 30˚, 45˚, 60˚, 75˚, and 90˚were performed using the RFPA 2D -Flow code.The list of rock and bedding parameters used in the simulations is presented in Table 3.The  As shown in Fig 6, the HF initiated and propagated with the greatest in situ stress at the beginning.When θ was 30˚(4σ = 2, 4, 6, 8 MPa), 45˚(4σ = 2, 4 MPa), and 60˚(4σ = 2 MPa), the HF dilated the BP and formed a simplex bi-wing fracture as the HFs approached the weak BP.When θ was 45˚(4σ = 6, 8 MPa), 60˚(4σ = 4, 6 MPa), 75˚(4σ = 2 MPa), and 90˚(4σ = 2 MPa), the BP shear slipped and formed a jagged fracture with a tortuous appearance.In this situation, a complex fracture network was formed by the frequent deflection of HFs and the penetration of the bedding.When θ was 60˚(4σ = 8 MPa), 75˚(4σ = 6, 8 MPa), and 90˚(4σ = 4, 6, 8 MPa), the HF penetrated the BPs and formed a horizontal HF, which propagated in the maximum principal stress direction.Overall, the smaller the intersection angle and the principal stress difference are, the easier the BP dilation by the HF is; the bigger the intersection angle and the greater the principal stress difference are, the easier the BP penetration by the HF will be.In addition, the weaker the bedding is, the easier the BP dilation by the HF will be.In general, the hydraulic fracture evolution was mainly controlled by the principal stress difference, the intersection angle, and the strength of the BP, which explains why these parameters are considered in the criterion.Young's modulus (E) GPa 30(meso value is 37) 3(meso value is 3. friction coefficient f = 1.21 due to the thickness of the BP [25].As shown in the figure, the dilation criterion and slip/penetration criterion curves were obtained after substitution.The simulations are consistent with the criterion.In Fig 7, the simulations and the criterion indicate that the HF would dilate the BP under low principal stress difference and intersection angle and penetrate the BP under high principal stress difference and intersection angle.The BP would shear slip under appropriate intersection angle and principal stress difference, which is conducive to the formation of a complex fracture network.

Numerical simulations for the influence of absolute value of principal stress on interaction
Considering the absolute value of principal stress would affect the stress state of the specimen, eight simulations with different maximum and minimum principal stress but same principal stress difference and intersection angle were performed using the RFPA 2D -Flow code.In this part,  As shown in Fig 8, as the maximum principal stress and the minimum principal stress increased from σ 1 = 4MPa and σ 3 = 2MPa to σ 1 = 18MPa and σ 3 = 16MPa uniformly, the fracture evolutions are basically unchanged.The fractures initiate and propagate in the maximum principal stress direction, and then the HF dilates the BP and forms a simplex bi-wing fracture as the HFs approach the weak BP.
It also can be seen from Fig 8 that the breakdown pressure of the specimens increases with the absolute values of the principal stress increasing although the fracture evolution does not change.In order to describe the relationship between the absolute values of the maximum and the minimum principal stress, a confining pressure ratio f s ¼ s 3 s 1 À s 3 is proposed.The larger f σ is, the smaller the confining pressure difference and the larger the absolute principal stress value become.As shown in Fig 9, the breakdown pressure of the specimens increases linearly with the increasing of the confining pressure ratio.In summary, it can be concluded that the absolute value of principal stress would not affect the interaction between BP and HF but influence the breakdown pressure of the specimen.

Discussions
In this study, a criterion is proposed to predict the interaction behaviors between HF and BP.In comparison with other criteria, the proposed prediction criterion considers three interactions (Dilate, Slip and Penetrate) but not only the penetration, and it is easy to use directly.In addition, other studies mostly focus on the interaction between HF and NF rather than HF and BP.Moreover, twenty simulations with different principal stress differences and intersection angles are simulated using the RFPA 2D -Flow code which is more suitable for rock than most of other finite element codes because it considers the heterogeneity of materials.The simulation results exhibit good agreement with the criterion results for a wide range of angles.In order to study the influence of the absolute value of principal stress on interaction, eight simulations with the same principal stress difference and intersection angle but different absolute values of the principal stress are simulated.The simulation results show that the absolute value of principal stress would not affect the interaction between BP and HF, which normalizes the application of prediction criterion.This study can provide relevant reference about the formation of complex fracture networks in shale layers.2. HFs propagate due to tensile stress, but the AE events showed that tensile failure is the main but not the only form of failure in hydraulic fracturing in the shale layer.Compression failure is rare but also exists.
3. The smaller the intersection angle and principal stress difference are, the easier the BP dilation by the HF will be; the larger the intersection angle and principal stress difference are, the easier the BP be penetration by the HF will be.The BP shear slip can form a complex fracture network through the deflection of HFs and the penetration of beddings, helping improve the reservoir permeability.

Fig 7
presents a comparison of the simulation results and the prediction criterion.According to the numerical model, the parameters in the criterion are taken as b a ¼ 3 (based on the simulation results in Fig 6, a = 0.04m and b = 0.12m), K IC ffi ffi ffi pa p ¼ 3:15 (K IC = 0.313+0.027Ewhich established by Whittaker in 1992 [36]), T O = 2 (rock's compressive-tensile strength ratio is taken as 10), S O = 2.5 (the shear strength is taken as 12.5% of the compressive strength), and